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Preface 


The  purpose  of  this  study  was  to  calculate  the  electric 
fields  in  the  air  resulting  from  a  surface  nuclear  burst. 

To  achieve  this,  a  computer  program  was  developed  which  used 
curve  fits  for  the  source  terms  and  the  electron  air  chem¬ 
istry  parameters  to  find  a  solution  to  Maxwell's  equations 
in  the  quasi-static  EMP  phase.  The  results  of  these  calcu¬ 
lations  were  compared  to  approximate  analytic  solutions. 

Since  the  solution  in  the  late-time  phase  is  time- 
independent,  the  computer  program  in  Appendix  F  should  be 
useful  because  of  the  relatively  fast  execution  time.  This 
allows  many  variations  in  parameters  to  be  studied  easily. 

I  would  first  like  to  thank  Maj .  Gary  Knutson  of  the 
Air  Force  Weapons  Laboratory  for  suggesting  this  project. 

A  great  deal  of  appreciation  must,  also,  be  extended  to 
my  faculty  advisor  Lt.  Col.  John  Erkkila.  He  guided  me  when 
things  went  poorly  and  gave  me  the  freedom  to  do  my  own 
work.  This  combination  greatly  enhanced  the  knowledge  I 
gained. 

Finally,  I  would  like  to  thank  my  wife,  Monica,  for  her 
patience  and  support  at  all  times.  This  thesis  truly 
belongs  to  both  of  us. 


James  R.  Downey 


Contents 


Page 

ii 


Preface  .... 

List  of  Figures .  v 

List  of  Tables  . viii 

Abstract . 

I.  Introduction  .  i 

Background  .  i 

Problem  .  3 

Scope  and  Assumptions .  3 

General  Approach  .  4 

Overview  .  4 

II.  Background  Quasi-Static  EMP  Theory  .  6 

Maxwell's  Equations  .  6 

Approximate  Solution  .  8 

Air  Chemistry .  10 

Limit  of  Approximate  Solution .  11 

III.  An  Improved  Analytic  and  Equivalent  Finite 

Difference  Solution  .  13 

Legendre  Polynomial  Expansion  .  13 

Numerical  Solution  Development  .  15 

Solution  Algorithm  .  17 

Boundary  Conditions  .  19 

Computer  Solution  .  20 

IV.  Initial  Numerical  Results  .  21 

Field  Calculations  .....  .  21 

Comparison  with  Expected  Conductivity  ....  23 

Comparison  with  Approximate  Expressions  ...  31 

Limits  of  Grover's  Model  .  33 

V.  Late-Time  EMP  Sources  and  Air  Chemistry 

Parameters .  34 

Late-Time  Sources  .  34 

Air  Chemistry  Parameters  .  41 


iii 


Page 


VI.  An  Improved  Numerical  Method  .  44 

Solution  Development  .  44 

Solution  Technique  .  51 

VII.  Results . 52 

Test  Problem  . . 52 

Theta  Dependence . 57 

Boundary  Condition  at  r  =  0 . 58 

Field  Dependence . 58 

Effect  of  JQ . 59 

Variation  in  Electron  Mobility  .  59 

Summary . 64 

VIII.  Parametric  Studies  .  67 

Code  Accuracy . 67 

Time  Dependence . 68 

Yield  Dependence . 70 

Water  Vapor  Dependence  .  74 

IX.  Conclusions  and  Recommendations  .  78 

Conclusions . 78 

Recommendations  .  79 

Bibliography  .  80 

Appendix  A:  Finite  Differences  in  Initial  Solution  .  .  82 

Appendix  B:  Initial  Numerical  Results  .  85 

Appendix  C:  Analytic  Fits  for  Sources  and  Air 

Chemistry  Parameters  .  94 

Appendix  D:  Romberg  Integration  .  98 

Appendix  E:  Parametric  Studies  .  101 

Appendix  F:  Program  Description  and  Listing  .  114 

Vita . I37 


1 


List  of  Figures 

Figure  Page 

1  Surface  Durst  Geometry  .  2 

2  Expected  Conductivity  Compared  to 

Grover's  Model  (3),  oq  =  178.0 .  25 

3  Expected  Conductivity  Compared  to 

Grover's  Model  (3),  =  143.2 .  26 

4  Er  vs  Radius,  Initial  Solution  .  27 

5  Ef  vs  Angle,  Initial  Solution  .  28 

6  -  Efl  vs  Radius,  Initial  Solution .  29 

7  -  EQ  vs  Angle,  Initial  Solution .  30  I 

8  Electronic  and  Ionic  Conductivity  vs  | 

Radius  for  Initial  Problem  .  32  f 


9  Ionization  Rate  vs  Radius  at  Four  Angles  ...  36 

10  Radial  Compton  Current  (-  Jr)  vs 


Radius  at  Four  Angles .  37 

11  Theta  Compton  Current  (J0)  vs 

Radius  at  Three  Angles  .  38 

12  -  J  and  JQ  vs  Time  .  40 

r  t) 

13  Electron  Attachment  Rate  vs  Total  Field  ...  42 

14  Electron  Mobil  t,  v*  Total  Field .  43 

15  Er  vs  Radius,  Full  Solution .  53 

16  Er  vs  Angle,  Full  Solution .  54 

17  -  E.  vs  Radius,  Full  Solution .  55 

18  -  E0  vs  Angle,  Full  Solution .  56 

19  E^  vs  Radius,  Field  Independent .  60 

20  -  E„  vs  Radius,  Field  Independent .  61 

21  E  vs  Radius,  J„  =  0 .  62 

r  o 


.v 


Figure 

22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 


List  of  Figures 
(continued) 

Page 


-  E0  vs  Radius,  Jy  =  0 .  63 

Er  vs  Radius,  Electron  Mobility  Not 

Limited .  65 

-  Ey  vs  Radius,  Electron  Mobility  Not 

Limited .  66 

Total  Field  vs  Time,  0  =  90° .  69 

Total  Field  vs  Time,  0  =  0° .  71 

Total  Field  vs  Yield,  0  =  90° .  72 

Total  Field  vs  Yield,  0  =  0° .  73 


Total  Field  vs  Water  Vapor  Content,  o  =  90°  .  75 
Total  Field  vs  Water  Vapor  Content,  0=0°..  76 


E  vs  Radius,  Simplified  Sources 

»o  =  178.0 .  86 

E  vs  Angle,  Simplified  Sources 

o  =  178.0 .  87 

o 

-  E0  vs  Radius,  Simplified  Sources 

o  =178.0 .  88 

o 

-  E()  vs  Angle,  Simplified  Sources 

oq  =  178.0 .  89 

E  vs  Radius,  Simplified  Sources 

=  143.2 .  90 

E  vs  Angle,  Simplified  Sources 

or  =  143.2 .  91 

o 

-  E A  vs  Radius,  Simplified  Sources 

o  =  143.2  . .  92 

o 

-  E„  vs  Angle,  Simplified  Sources 

o  0  =  143.2 .  93 

Generalized  Romberg  Terms  .  99 

Romberg  Algorithm  .  100 


vi 


•v 


List  of  Figures 
(continued) 


1 

Figure 

Page 

! 

41 

Total  Field  vs 
t  =  10~2  sec  . 

Radius,  Pull  Solution 

.  .  102 

i 

i 

| 

42 

Total  Field  vs 
t  =  10~2  sec  . 

Angle,  Full  Solution 

•  •  103 

f 

43 

Total  Field  vs 
t  =  10-1  sec  . 

Radius,  Full  Solution 

•  •  104 

44 

Total  Field  vs 
t  =  10-1  sec  . 

Angle,  Full  Solution 

•  •  105 

| 

45 

Total  Field  vs 
Y Id  =  1000  KT 

Radius,  Full  Solution 

.  .  106 

i 

46 

Total  Field  vs 
Yld  =  1000  KT 

Angle,  Full  Solution 

•  •  ]  07 

! 

47 

Total  Field  vs 
Yld  =  100  KT  . 

Radius,  Full  Solution 

•  •  108 

48 

Total  Field  vs 
Yld  =  100  KT  . 

Angle,  Full  Solution 

•  •  109 

( 

49 

Total  Field  vs 
Wvc  =0.00  .  . 

Radius,  Full  Solution 

•  •  110 

t 

b 

« 

50 

Total  Field  vs 
Wvc  =0.00  .  . 

Angle,  Full  Solution 

•  •  111 

51 

Total  Field  vs 
Wvc  =0.04  .  . 

Radius,  Full  Solution 

.  •  •  112 

52 

Total  Field  vs 
Wvc  =0.04  .  . 

Angle,  Full  Solution 

.  .  .  113 

t- 

53 

Program  Flow  Diagram  . 

.  .  .  115 

vn 


List  of  Tables 


Table  page 

I  Computer  Code  Parameters  Used  in 

Calculations . . 68 

II  Analytic  Source  Functions  .  95 

III  Mobility  and  Attachment  Rate  Equations  ....  97 


viii 


AFIT/GNE/PH/83M-5 


Abstract 

A  numerical  solution  was  developed  to  find  the  quasi¬ 
static  EMP  fields  resulting  from  a  surface  nuclear  burst. 

By  ignoring  the  time  derivatives  in  Maxwell's  equations  and 
expanding  the  electrostatic  potential  in  Legendre  polyno¬ 
mials,  a  set  of  differential  equations  was  obtained  that  was 
dependent  on  r  only  (r  =  radius  from  burst).  By  mploying 
finite  differences,  a  tri-diagonal  matrix  equate  was 
obtained  that  was  non-linear  in  the  electric  fie?  This 
equation  was  solved  using  an  iterative  scheme. 

Analytic  curve  fits  based  on  Monte  Carlo  data  were  used 
to  determine  the  air  ionization  rate  and  Compton  current. 
This  permitted  the  angular  variation  of  these  sources  to  be 
included  in  the  calculations.  In  addition,  the  electron  air 
chemistry  parameters  were  allowed  to  vary  with  the  local 
electric  field  and  water  vapor  content  of  the  air.  This  was 
accomplished  by  the  use  of  previously  developed  equations 
for  these  terms. 

Field  calculations  were  made  and  compared  to  previous 
analytic  results.  The  differences  were  found  to  be  as  great 
as  50  percent  for  peak  field  values.  A  limited  parametric 
study  was  performed  to  consider  the  effects  of  varying  the 
time,  the  yield,  and  the  water  vapor  content.  The  computer 
program  documented  in  this  report  should  be  useful  for  late- 


time  EMP  calculations  because  of  the  short  execution  time. 
Furthermore,  the  code  was  written  such  that  improvements  in 
the  various  fits  could  easily  be  incorporated. 


THE  CALCULATION  OF  LATE-TIME 


SURFACE  BURST  EMP  FIELDS  USING 
A  TIME- INDEPENDENT  NUMERICAL  METHOD 

I.  Introduction 


Background 

In  recent  years  interest  has  increased  in  determining 
the  electromagnetic  pulse  (EMP)  produced  by  surface  nuclear 
bursts.  Of  special  concern  is  the  EMP  that  occurs  in  the 
late-time  regime  (-  10~‘'  to  1  seconds,  retarded  time).  This 
has  largely  been  motivated  by  changing  system  requirements 
and  a  desire  to  explain  such  things  as  the  lightning  strokes 
observed  during  the  MIKE  test  shot  (Ref  16). 

When  a  surface  nuclear  burst  occurs,  electric  currents 
are  set  up  in  the  air  as  a  result  of  Compton  scattered 
electrons.  The  gamma  rays  penetrating  the  ground  are  highly 
attenuated;  thus,  the  EMP  source  region  is  basically  a 
hemisphere  above  the  surface.  The  generated  Compton  cur¬ 
rents  are  primarily  in  the  radial  direction  producing  a 
radial  electric  field;  however,  the  presence  of  the  ground 
creates  transverse  field  components  which  can  exist  well 
beyond  the  source  region  (Ref  6:31-32).  The  geometry  for 
the  source  region  is  shown  in  Figure  1. 

The  problem  of  late-timc  EMP  has  been  investigated  by 
several  authors  (Ref  4;  8;  9;  and  17).  The  basis  of  the 
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Figure  1.  Surface  Burst  Geometry 

theory  assumes  that  the  time  derivatives  in  Maxwell's 
equations  will  he  small  at  times  greater  than  approximately 
100  us  and  can  be  ignored.  This  assumption  is  known  as  the 
quasi-static  approximation.  The  electric  field  is  then 
derivable  from  a  scalar  potential  F.  =  -qrad  . 

Longmire  (Ref  8),  Hill  (Ref  4),  and  Wyatt  (Ref  17) 
assumed  that  the  radial  component  of  the  electric  field 
would  be  much  smaller  than  the  polar  component;  therefore 
they  ignored  it  in  their  calculations.  More  recently, 
Grover  (Ref  2)  considered  the  problem  of  quasi-static  EMP 
without  making  any  assumptions  about  the  field  components. 
Using  simplified  air  conductivity  models  and  expanding  the 
electric  potential  in  Legendre  polynomials,  Grover  obtained 
analytic  solutions  for  the  electric  fields. 

Though  Grover's  analytic  results  are  useful,  the  solu¬ 
tions  are  limited  by  several  factors:  (1)  the  simplified 
conductivity  models  Grover  used  are  only  valid  for  certain 
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space-time  regimes  (Ref  2:23-25),  (2)  the  source  terms  (air 
ionization  rate  and  radial  Compton  current)  are  restricted 
because  no  polar  variation  is  allowed,  and  (3)  the  electron 
mobility  and  attachment  rate  are  assumed  to  be  independent 
of  the  electric  field  and  water  vapor  content  of  the  air. 
These  assumptions  are  necessary  to  obtain  analytic  results. 

Problem 

The  problem  investigated  in  this  study  is  to  determine 
the  magnitude  of  the  electric  field  in  the  air  during  the 
quasi-static  phase  resulting  from  a  surface  nuclear  burst. 

The  angular  dependence  of  the  source  terms  and  varia¬ 
tions  in  the  air  chemistry  parameters  are  included  in  the 
calculations . 

Scope  and  Assumptions 

The  time  range  of  interest  in  this  report  is  from  10-4 
to  10  1  seconds.  This  range  includes  the  time  frame  of  the 
reported  nuclear  lightning  (=  1  ms)  and  allows  for  the 
dominance  of  ground  or  air  capture  sources. 

The  conductivity  of  the  ground  is  usually  much  larger 
than  that  of  the  air,  especially  for  later  times  and  smaller 
yields;  therefore  it  will  be  considered  infinite.  This  will 
permit  the  use  of  simplified  boundary  conditions  for  the 
solution  of  Maxwell's  equations. 

Self-consistent  effects  between  the  generated  fields 
and  the  source  currents  are  ignored  in  all  calculations. 

This  is  done  to  simplify  the  solution  technique. 
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General  Approach 

The  first  phase  in  this  study  was  the  development  of  a 
numerical  solution  which  did  not  include  the  polar  variation 
of  the  sources  or  any  variations  in  the  air  chemistry  param¬ 
eters.  This  allowed  a  comparison  and  confirmation  of  the 
results  presented  by  Grover  and  provided  a  basis  for  com¬ 
parison  with  later  calculations .  The  next  step  was  to 
derive  an  angular  and  field  dependent  solution  (ignoring 
self-consistency)  using  curve  fits  for  the  source  terms 
(Ref  10)  and  air  chemistry  parameters  (Ref  7) .  The  results 
obtained  here  were  then  compared  to  the  earlier  calcula¬ 
tions.  Finally,  limited  parametric  studies  were  performed 
to  consider  the  effects  of  varying  the  time,  the  weapon 
yield,  and  the  water  vapor  content. 

Overview 

In  Chapter  II  the  background  theory  of  quasi-static  EMP 
is  discussed.  The  development  of  the  initial  solution  using 
the  simplified  models  is  presented  in  Chapter  III.  The 
results  of  calculations  using  the  first  solution  technique 
and  a  comparison  with  one  of  Grover's  models  are  given  in 
Chapter  IV.  In  Chapter  V  the  angular  dependence  of  the 
source  terms  and  the  field  and  water  vapor  dependence  of  the 
air  chemistry  parameters  are  discussed.  The  field  and  polar 
dependent  derivation  and  solution  algorithm  are  given  in 
Chapter  VI.  Chapter  VII  presents  a  comparison  of  the 
results  obtained  using  the  new  model  with  the  previous 
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calculations.  A  series  of  parametric  studies  are  given  in 
Chapter  VIII  and  conclusions  and  recommendations  are 
presented  in  Chapter  IX. 
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II.  Background  Quasi-Static  EMP  Theory 

By  making  a  number  of  simplifying  approximations,  the 
quasi-static  EMP  fields  can  be  found  by  analytic  solutions 
to  Maxwell’s  equations.  In  this  chapter  the  basic  theory 
behind  these  approximate  results  is  presented.  The 
limitations  of  this  theory,  also,  are  discussed. 


Maxwell ' s  Equations 

As  mentioned  in  the  introduction,  it  is  commonly 
assumed  that  at  times  greater  than  approximately  100  ps 
after  a  surface  nuclear  burst  the  electric  fields  become 
quasi-static.  Consequently,  the  time  derivatives  in 
Maxwell's  equations  are  small,  compared  to  the  other  terms, 
and  can  be  set  equal  to  zero  (Ref  9:44).  Hence, 

V  x  £  =  0  (2.1) 

_L_ 

w Q  $  x  S  =  +  3C  (2.2) 

where  2  is  the  electric  field  (volts/m) ,  §  is  the  magnetic 
field  (webers/m; ) ,  o  is  the  conductivity  of  the  air  (mho/m), 
Jc  is  the  source  current  density  (Amps/m2),  and  yQ  is  the 
magnetic  permeability  (henries/m) . 

From  vector  analysis  it  is  known  that  $  x  (v$)  =0, 
where  41  is  some  scalar;  therefore,  E  can  be  defined  as 

E  =  -v*  (2.3) 
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where  <t>  is  the  scalar  electric  potential  (volts)  .  The 
individual  field  components  are  found  by 


Er  =  - 


341 

(2.4) 

3r 

a  <t> 
ae 

(2.5) 

where  r  is  the  radius  from  the  burst  (meters) .  E,  is  zero 

(h 

i 

from  symmetry. 

The  permeability  in  Eq  (2.2)  is  assumed  to  be  constant. 
Thus,  the  equation  can  be  rewritten  such  that 

V  •  (If  xS)  =0  (2.6) 

Mo 

or 

$  *  +  Jc)  =  0  (2.7) 

Substituting  Eq  (2.3)  into  Eq  (2.7)  gives 

V  .  ( u  7  4> )  =  V  •  J  (2.8) 

c 

Expanding  Eq  (2.8)  in  spherical  coordinates  and  assuming 
azimuthal  symmetry  leads  to 


-L  -i-  (r2  ii,  +  _L  _J_  _L  (SING  i±)  +  I  ii 
r2  or  3r  r2  SINO  30  30  a  3r  3r 

(2.9) 

+  -L-  ii  ii  =  _i — L  ( r 2  j  )  +  — I —  -i.  (siNe  j  ) 

or2  36  30  or2  3r  C  or  SINO  36 

This  is  the  equation  that  must  be  solved  to  find  the 
electric  field  in  the  air. 


Approximate  Solution 


Longmire  (Ref  9) ,  Hill  (Ref  4) ,  and  Wyatt  (Ref  17)  have 

shown  that  an  analytic  solution  to  Eq  (2.9)  can  be  found  if 

some  simplifying  assumptions  are  made. 

Since  E  * <  E  at  late  times  it  is  set  equal  to  zero, 
r  a 

In  addition  the  air  conductivity  and  radial  Compton  current 
are  assumed  to  be  independent  of  the  polar  angle  theta. 
Finally,  Monte  Carlo  studies  have  shown  (Refs  5  and  10)  that 
the  theta  component  of  the  Compton  current  source  is  much 
smaller  than  the  radial  component,  except  for  small  r;  thus 
J()  is  ignored. 

Using  all  of  the  above  assumptions  Eq  (2.9)  reduces  to 

- 1 - L  (si No  H  )  -  — - -  ( r •’  J  )  (2.10) 

r  SINO  96  or2  9r  r 


Longmire  has  shown  (Ref  8)  that  the  radial  Compton 
current  can  bo  fairly  well  approximated  by 


-r/\ 

J  =  J  -  (Amps/m2) 

r  o  o 


(2.11) 


where  Jq  is  a  constant  for  a  given  time  and  yield,  and  X  is 
the  effective  gamma  ray  mean  free  path  (meters) . 

Substituting  this  equation  into  Eq  (2.10)  leads  to 


— - —  (SINO  k  )  =  —  J  e"r/X 

SINO  30  oX 


(2.12) 


or 

— - —  (SINO  E  )  =  —  J  (2.13) 

SINO  30  0  oX 


•v 
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Integrating  with  respect  to  0  gives 


SINO  E. 

O 


-  —  (l  -  cose) 

A  a 


(2.14) 


or 

E.  =  —  — -  TAN  -  (2.15) 

6  A  a  2 

This  is  Longmire's  result.  By  including  JQ ,  Wyatt  obtained 
the  following  (Ref  17:8) 

J 

E  =  -  (  — £  TAN  -  -  JA)  (2.16) 

u  o  o 

0  A  2. 

which  reduces  to  Eq  (2.15)  if  J  =  0. 

0 

From  Eq  (2.1)  and  the  assumption  that  =  0 ,  one  sees 
that  E  must  satisfy  the  following 

-  -  (r  E  )  =0  (2.17) 

r  3r 

or 

E  <*  —  (2.18) 

9  r 

For  the  above  condition  to  be  true  the  ratio  of  Jr  over 
o  must  be 


(2.19) 


As  will  be  shown  below,  this  requirement  places  a 
restriction  on  the  range  over  which  Eq  (2.15)  applies. 


Air  Chemistry 

To  solve  for  the  electric  field,  the  conductivity  must 
be  known.  The  total  conductivity  consists  of  contributions 
from  free  electrons  and  positive  and  negative  ions. 

The  electron  conductivity  can  be  written  (Ref  6:21-22) 

o  =  e  ii  —  (mho/m)  (2.20) 

e  e  a 


where  e  is  the  charge  on  an  electron  (coulombs) ,  pe  is  the 
electron  mobility  (m7/V*sec),  S  is  the  local  ionization  rate 
(ion  pairs/m3 -sec) ,  and  is  the  electron  attachment  rate 
(sec-1).  Both  mg  and  are  dependent  on  the  total  field, 
but  for  now  are  considered  field  independent. 

At  late  times  most  of  the  electrons  are  attached;  thus, 
the  number  of  positive  ions  (N+)  is  approximately  the  same 
as  the  number  of  negative  ions  (N_)  and  the  ionic  conduc¬ 
tivity  can  be  written  as 


o 


I 


2 


(mho/m) 


(2.21) 


where  m  r  is  the  ion  mobility  (m''/V*sec)  and  is  the  ion- 
ion  recombination  rate  (m3/sec) .  The  two  represents  the 
fact  that  N+  =  N_  . 

The  ionization  rate  in  Eqs  (2.20)  and  (2.21)  is  assumed 
(Ref  8)  to  be  of  the  form 


S  = 


(2.22) 
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where  Sq  is  a  constant  for  a  given  time  and  yield. 

The  total  conductivity  can  then  be  written  (using  Eq 
(2.22)  for  S) 


(2.23) 


As  seen  in  the  radial  dependence,  electrons  should 
dominate  the  conductivity  close  to  the  burst  and  ions 
farther  out. 


Limit  of  Approx imate  Solution 

If  electrons  dominate  the  conductivity,  J  /a  is 
nearly  independent  of  r.  Thus  it  would  be  expected  that 
the  Longmire  model  would  be  inaccurate  at  shorter  ranges 
from  the  burst.  If  ions  dominate,  the  ratio  becomes 


r 

a 


-r/2  A 


(2.24) 


or 

re  a  _  constant  (2.25) 

Longmire  (Ref  9:57)  states  that  this  requirement  is  not 
accurately  valid  and  he  limits  the  range  of  applicability  to 
0.2  to  2.0  kilometers  with  a  maximum  variation  of  +  30 
percent . 


EMP  theory  predicts,  and  others  have  shown  (Ref  13), 
that  the  electron  conductivity  will  be  important  even  out  to 


approximately  1000  meters.  This  obviously  presents  a 
problem  for  the  model  discussed  above.  In  the  next  chapter 
an  analytic  solution  will  bo  derived  which  is  not  limited 
by  the  above  requirements. 
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III.  An  Improved  Analytic  and 
Equivalent  Finite  Pi f ference  Solution 

In  the  last  chapter  an  analytic  solution  to  Maxwell's 
equations  was  found  by  making  several  simplifying  assump¬ 
tions.  Recently,  Grover  (Ref  2)  obtained  an  improved 
solution  by  expanding  the  electric  potential  in  Legendre 
polynomials  and  by  using  simplified  models  for  the  air 
conduct i vi ty . 

In  this  chapter  the  basis  of  Grover's  analytic  solution 
is  presented  along  with  an  equivalent  finite  difference 
method.  The  numerical  solution  was  developed  to  provide  a 
basis  for  comparison  with  the  analytic  results  (Ref  2)  and 
the  more  general  results  presented  later  in  this  report. 

Legendre  Polynomial  Expansion 

Longmire  and  Wyatt  assumed  that  the  radial  component  of 
the  field  would  be  smaller  than  the  theta  component  and  they 
ignored  it  in  their  calculations .  Grover  (Ref  2),  however, 
found  an  analytic  solution  without  making  such  an  approxi¬ 
mation  . 

By  expanding  the  potential  in  a  series  of  Legendre 
polynomials,  a  set  of  ordinary  differential  equations  was 
obtained  which  only  depended  on  r  (Ref  2:7).  The  potential 
is  written  as 

=  Z  A  (r)  P  (x) 
t=odd  l 
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<|>  ( r ,  o ) 


(3.1) 


where  A  (r)  are  unknown  coefficients  to  be  determined,  and 
I*  ( x )  are  the  Legendre  polynomials  (x  =  COSu).  The  angle 
theta  is  measured  from  the  vertical  such  that  u  =  0°  points 
straight  up  and  >>  =  90°  lies  along  the  ground. 

Choosing  the  summation  over  odd  it's  satisfies  the 
boundary  conditions  on  the  field:  namely,  that  the  radial 
component  of  the  electric  field  must  be  zero  at  0  =  90°  for 
an  infinite  conductivity  ground  plane,  and  the  polar  com¬ 
ponent  must  be  zero  at  o  =  0°  due  to  symmetry. 

To  solve  the  problem  analytically,  Grover  assumed  that 
the  radial  Compton  current  and  the  ionization  rate  were 
independent  of  theta.  The  expressions  presented  earlier, 

Eqs  (2.11)  and  (2.22),  were  used  for  these  parameters. 

Instead  of  using  the  actual  air  chemistry,  Grover 
approximated  the  conductivity  by  three  simplified  models. 
Model  (1)  was  designed  to  represent  the  early-time  regime 
when  the  conductivity  is  dominated  by  electrons.  Model  (2) 
approximated  the  late-time  regime  when  ions  dominate.  Model 
(3)  was  used  for  intermediate  times  when  electrons  dominate 
in  close  and  ions  farther  out  (Ref  2:7).  These  models  were 
necessary  to  obtain  analytic  solutions. 

The  purpose  of  this  work  was  not  to  get  Grover’s 
analytic  solution;  consequently  their  derivation  is  left  to 
the  reader  (Ref  2) .  Rather,  an  equivalent  numerical  method 
is  presented  below.  This  solution  was  derived  to  provide  a 
basis  for  comparison  to  Grover’s  results  and  act  as  a 
stepping  stone  to  the  more  general  results  presented  later. 
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Numeric,!  1  Solut  ion  Devc  lopment 

Recalling  the  differential  equation  presented  in  Eq 
(2.9)  and  substituting  Eq  (3.1)  for  the  potential  leads  to 


-odd 


'•A  (r) 

P  (x) - - -  + 

3r2 


I2-  +  i 

\r  o  Jr/ 


E  P.(x) 

c=odd 


3A£(r) 


3r 


(3.2) 


-  —  (  SINu 

r'SINO  .)•>  \ 


«-  \  i  a  J  r  2 

>:  A  (r)  P  (x)  )  =  -  — -  +  —  J 
i  =odd  /  o  3r  or 


where  J  is  assumed  to  be  zero. 

1 

From  Legendre's  differential  equation  it  can  be  shown 
(Ref  14:134)  that 


C_*  I  M  i  .  ' 


;»p  (x) 


SINO 


S  I  N 


I  (I 


)  =  -t(Hl)  Pa(x)  (3.3) 


Substituting  this  relationship  into  Eq  (3.2)  and  rearranging 
the  terms  gives 


):  P  (x) 

e=odd 


'  A‘(rl  +| 

f  2  1  3  o  ^ 

3Ae(r)  lU+llA^r)  ' 

3r  2 

^  r  o  3r/ 

3r  r2 

(3.4) 


-  - -  +  —  J 

o  3r  or  1 


From  the  orthogonality  of  Legendre  polynomials,  the 
following  is  obtained  (Ref  14:131) 


+  1 

/  P  (x)  P  (x)  dx 
J  m  n 

-1 


0  m  f  n 

2 

~ — pr  m  =  n 
2n+l 


(3.5) 
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For  odd  or  even  polynomials 


1  +1  0  m  ^  n 

fp _(x)  P  (x)  dx  =  i  /Pm<x)  p  (x)  dx  =  (3.6) 

J  m  n  2  ./  m  n 

0  -1  1 


2n+l 


Therefore,  the  summation  in  Eq  (3.4)  can  be  eliminated 

by  multiplying  both  sides  of  the  differential  equation  by 

I'  (x)  and  integrating  from  0  to  1 .  Thus 
m 


-  Him  a  ,n 

;-r'  \r  a  :>r  /  3r  rJ 


(3.7) 


=  (2  P  +  1 )  J  P  (x)  f  -  — -  +  —  J  ]  dx 


o  'J  r  or 


The  term  in  brackets  on  the  right-hand  side  of  the 
equation  can  be  removed  from  the  integral  because  Jr  and  o 
do  not  depend  on  theta.  Consequently,  Eq  (3.7)  becomes 


3  A  j  (  r ) 


.  +  1  i!  \  ^«lrl  -  «nn>  a  (r , 

\r  o  >r  J  jr  r? 

,  2j  i  y 

+  1)  1  — 1  +  — -  \  J  P  e  (x)  dx 
.  o  3r  or  J  0 


(3.8) 


Note  that  the  angular  dependence  which  originally 
existed  in  the  differential  equation,  Eq  (2.9),  has  been 
integrated  out  and  a  series  of  one-dimensional  equations 
remain.  Though  different  from  Grover's  radial  equations, 
the  solution  should  yield  the  same  results.  Since  Eq  (3.8) 
is  solved  numerically,  it  has  the  added  advantage  of  not 


requiring  a  simplified  model  for  the  conductivity.  The 
assumed  air  chemistry  equations  can  be  used. 


Solution  Algorithm 

To  solve  Eq  (3.8),  finite  differences  were  used. 

Letting  r^  =  n-Ar  where  n  is  the  n ' th  radial  position  where 
fields  are  to  be  found,  all  the  derivatives  in  the  equation 
are  replaced  by  their  equivalent  central  difference  approxi¬ 
mations.  It  can  be  shown  (see  Appendix  A),  after  substitu¬ 
tion  of  the  difference  operators  and  a  rearrangement  of 
terms,  Eq  (3.8)  becomes 


a  A  .  +  b  A  +  c  A„  .  =  R„  (3.9) 

n  ?,n+l  v,n  £,n  n  £,n-l  £,n 


where  a  ,  b  ,  c  ,  and  R  are  terms  that  can  be  determined 
v.  / 1 1  n  x.  j  it 

at  any  range  point  n.  Hence,  for  each  given  £  the  equation 

is  solved  to  find  A  at  all  range  points. 

f  n 

Eq  (3.9)  is  a  tri-diagonal  matrix  equation.  The  method 

that  was  used  to  solve  it  was  a  two-sweep  algorithm  sometimes 

known  as  the  Thomas  method.  The  algorithm  is  outlined  below. 

Assume  that  the  coefficients  A  can  be  written  such 

£ ,  n 

that  A „  ,  ,  is  a  function  of  the  previous  coefficient  A„  , 

i .  e . 


Likewise 


A 


£  ,n  +  l 


£ ,  n 


e  + 
n 


f 

n 


(3.10) 


A.  =  A.  .  •  e  ,+f  , 

£,n  £,n-l  n-1  n-1 

Substituting  Eq  (3.10)  into  Eq  (3.9)  gives 
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(3.11) 


* 


s* 
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Approx ima t i nq 


by  its  finite  difference  equivalent 


(r) 


jives 


.)r 


E  =  -  >'• 

rn  .  =odd 


,n+l  A£,n-1 

2  Ar 


Vx) 


(3.18) 


Similarly,  E  is  found  by 


Eo  =  - 


ilk,  l  ?  A  (r)  P  *  (x) 
;)o  r  r  £=odd  £ 


(3.19) 


or 


E“  =  ~  L,  A«.n  V*1 

n  r  '’.=odd  ' 

n 


(3.20) 


where  P,(x)  is  the  l'st  Associated  Legendre  polynomial. 

Once  again  x  =  COS(O).  By  symmetry  the  azimuthal  component 

( E  )  is  zero. 

■t1 


Boundary  Cond i t ions 

To  find  all  the  f  's  and  e  ' s  in  Eqs  (3.15)  and  (3.16), 

n  n 

boundary  conditions  must  be  applied  on  the  fields  at  r  =  0 
and  r  =  r^  (N  =  maximum  number  of  range  points) . 

At  r  =  0  the  electric  field  must  be  finite.  Inspection 
of  Eq  (3.19)  shows  that  this  boundary  condition  is  satisfied 
if 


<|>  =  0  at  r  =  0 


(3.21) 


Thus,  from  Eq  (3.10) 


A 


£,1 


e  +  f 


o  o 


(3.22) 
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t 


or 


A 


2,1 


f 

o 


(3.23) 


As  r  increases  to  its  maximum  value  it  is  assumed  that 
the  radial  field  approaches  zero  (actually  Er  =  0  at  r  =  « ) . 
From  Eq  (3.17)  one  sees  that  this  condition  is  met  if 

*N+1  =  *N  (3*24) 


Hence 


asl,n+1  ah,n 


=  A 


2,N 


'N 


+  f. 


(3.25) 


which  requires  f^  =  0  and  e^  =  1 . 

With  the  boundary  conditions  satisfied,  a  complete 
solution  is  obtained. 


Computer  Solution 

Applying  the  derivation  and  the  algorithm  presented 
above,  a  program  was  written  and  used  to  find  the  electric 
field  assuming  simplified  sources.  The  results  of  these 
initial  numerical  calculations  are  presented  in  the  next 
chapter . 


r 
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IV.  Initial  Numerical  Results 


Using  the  alqorithm  developed  in  the  last  chapter  and 
simplified  sources,  the  electric  fields  were  calculated  for 
a  specific  yield/time  combination.  The  results  are  discussed 
in  this  chapter  and  compared  to  Grover's  analytic  solution 
for  the  same  problem.  In  addition,  a  comparison  is  made 
between  the  solution  using  Grover's  Model  (3)  conductivity 
(Ref  2:13)  and  the  approximate  air  chemistry  equations  pre¬ 
sented  earlier. 

r  i  <  •  Id  (.'a  L  eu  1  a  L  i  ons 

To  solve  for  the  electric  field,  the  source  terms  and 
air  chemistry  parameters  must  be  known.  The  particular  case 
investigated  was  a  10  megaton  surface  burst.  The  field  was 
calculated  at  10"  1  sec  (retarded  time)  after  the  burst.  This 
example  was  chosen  to  allow  a  direct  comparison  to  Grover's 
calculations  (Ref  2:25). 

The  ionization  rate  is  found  from  Eq  (2.22),  i.e. 


A  value  of  SQ  =  1.1  x 103°  ion-pairs/m *sec  was  chosen 
for  this  problem  (Ref  2:25).  320  meters  was  used  for  A 

(Ref  2:20) . 

The  radial  Compton  current  is  approximately  proportional 
to  the  local  ionization  rate  (Ref  8) .  Thus,  Jr  can  be 
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written  as 


J  =  c  •  s 
r 


(4.2) 


where  C  is  the  constant  of  proportionality.  The  value  of  C 
used  was  C  =  -8.2  x  lO"''1*.  This  number  was  based  on  an 
average  of  various  reported  values  (Ref  2:20).  Recalling 
Eq  (2.11)  for  Jf/  it  can  be  seen  that 


J  =  C 
o 


(4.3) 


Hence,  lor  the  example  considered  here 


J  =  -9.02  x  10 
r 


,-r/A 


(Amps/m2 ) 


(4.4) 


As  discussed  earlier,  Grover  used  simplified  models  for 
the  conductivity  instead  of  the  approximate  air  chemistry 
equations.  In  the  calculations  presented  in  this  work, 
Grover's  Model  (3)  was  used  (Ref  2:13).  This  model  was 
designed  to  provide  the  best  fit  for  the  conductivity.  Using 
this  model,  o  is  found  as 


-r/2  A 


(mho/m) 


(4.5) 


where  o  is  determined  by  a  fit  to  the  expected  true  form 


o-o  +  o  T  . 
e  I 


In  his  calculations,  Grover  assumed  (Ref  2:25): 


—  =  -63  kV/m 
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Thus,  from  I'.q  (4.4)  oQ  =  143.2.  It  was  found,  however, 

that  this  number  was  inconsistent  with  a  plot  of  the 

function  shown  in  Grover's  report  (Ref  2:14).  From  the 

curve  an  average  value  of  o()  =  178.0  was  estimated.  The 

discrepancy  could  not  be  explained;  therefore  calculations 

were  made  usinq  both  values  for  a  . 

o 

Plots  of  the  radial  and  polar  electric  fields  that  were 

found  usinq  the  above  sources  and  o  =  178.0  are  qiven  in 

o 

Appendix  B.  The  answers  aqreed  quite  well  with  those  found 
by  Grover;  varyinq  only  a  percent  or  two  for  shorter  ranges. 
The  difference  is  most  likely  due  to  the  uncertainty  in  the 
v,i  hie  <jI  -1  .  C  ileulat  ions,  also,  wore  made  usinq  a  =  143.2. 
These  results  are  given  in  Appendix  B.  The  fields  found 
using  this  value  of  oQ  were  slightly  larger  than  the  ones 
obtained  by  Grover. 

Comparison  wi th  Expected  Conductivity 

The  problem  encountered  in  the  above  calculations 
points  out  one  of  the  limiting  factors  of  the  simplified 
conductivity  models;  the  value  of  must  be  determined  from 
the  expected  true  form  of  the  conductivity.  From  Eqs  (2.21) 
and  (2.22)  it  was  shown  that 


As  will  be  seen  in  the  next  section,  the  electron 
mobility  and  attachment  rate  are  dependent  on  the  total 
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electric  field  and  water  vapor  content  of  the  air.  However, 
the  ratio  Me/ac  is  not  as  dependent  for  lower  water  vapor 
contents  and  these  parameters  are  often  approximated  by 
average  values. 

Shown  in  Figures  2  and  3  are  plots  of  Grover's  Model 
(3)  conductivity  (using  aQ  =  178.0  and  oq  =  143.2)  and  the 
expected  true  form  given  by  Eq  (4.6).  The  following  average 
values  were  assumed  for  the  air  chemistry  parameters  (Ref  2; 
8;  9) 

=  0.25  (m2/V*sec) 

a  =  1.5  x  108  (sec-1) 

e 

=  2.5  x  10-4  (m2/V-sec) 

Yr  =  2.0  x  10~12  (m3/sec) 

These  values  have  a  reported  uncertainty  of  ±  30  to  50 
percent  ( Re f  2)  . 

As  can  be  seen  in  both  figures,  the  model  provides  a 
good  estimate  of  the  expected  conductivity  over  much  of  the 
range . 

Since  the  calculations  made  in  this  report  were  done 
numerically,  Eq  (4.6)  could  be  used  to  find  the  conductivity 
instead  of  Grover's  model.  The  results  of  these  calcula¬ 
tions  are  shown  in  Figures  4-7.  The  fields  found  here  were 
slightly  larger  than  those  obtained  by  Grover.  This  might 
be  expected  since  the  conductivity  is  lower  in  the  middle 
part  of  the  range  (Figures  2-3)  .  The  model  appears,  however. 


Conductivity  (mho/m) 


Fiyure  3.  Expected  Conductivity  Compared  to 
Grover’s  Model  (3),  oq  =  143.2 
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LEGEND 


Figure  4.  E  vs  Radius,  Initial  Solution 
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Figure  5.  E  vs  Angle,  Initial  Solution 


15.0  30.0  450  600  750  900 

Angle  (deg) 

Figure  7.  -  EQ  vs  Angle,  Initial  Solution 

Eg.  2.15  Shown  as  a  Dashed  Line  (r = 1500m) 


to  predict  quite  well  the  electric  fields  that  would  be 
expected  if  the  air  chemistry  equations  were  used. 


Conipar ison  with  Approximate  Expressions 

From  the  work  of  Longmire  and  Wyatt,  it  was  found  that 
if  Ef  was  ignored,  an  analytic  expression,  Eq  (2.15),  for 
Eq  could  be  derived.  This  approximate  result  is  compared 
to  the  Grover-type  solution  in  Figures  6  and  7.  The  major 
difference  between  the  two  derivations  was  the  inclusion  of 
E^  by  Grover.  He  concluded  (Ref  2) : 

"...  it  is  found  that  the  neglect  of  the 
radial  component  of  the  electric  field  is  not 
justified  when  electrons  contribute  significantly 
to  the  air  conductivity  ...  even  in  the  special 
case  where  ions  dominate  the  air  conductivity,  the 
previous  approximate  solutions  . . .  are  only 
qualitatively  correct."  (Ref  2:7-8) 

In  Figure  8  the  electronic  and  ionic  components  of  the 
conductivity  from  Eq  (4.6)  are  plotted  as  a  function  of 
radius  from  the  burst.  It  is  seen  that  electrons  do  indeed 
dominate  siqma  out  to  approximately  1000  m;  hence  their 
effect  can  not  be  ignored. 

As  shown  in  Chapter  II,  E  was  required  to  be  <*  to  1/r. 
This  came  from  the  equation  for  the  curl  of  2,  i.e. 


1 

r 


3  r 


r  E, 


0 


(4.7) 


where  in  the  previous  work  ^E^./ 38  was  assumed  small,  and 
consequently,  ignored.  This  proved  to  be  incorrect  if 
electron  conductivity  was  important  and  partially  limited 
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-  it* 


Radius  (m) 


Fiqure  8.  Electronic  and  Ionic  Conductivity  vs 
Radius  for  Initial  Problem 
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for  ion  dominated  conductivity.  By  contrast,  Eq  (4.7)  is 
satisfied  exactly  by  Grover's  model.  Recalling  the  defining 


equations  for  H 

and  l 

4j  t 

0 

Eqs 

(3.17)  and 

(3.19)  , 

one  sees 

i 

-( 

r 

- 1  D«J> 

)  +  JL  M.  ' 

=  0 

(4.8) 

r 

3i  ' 

r  DO 

7  36  3r 

or 

3 

3 

3  3,, 

(<J>) 

- - -  (♦> 
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In  other  words  the  term  involving  E  is  the  same  as  the 

r 

one  with  . 

L 1  in  i_Ls  oj  Grover  '  s  Model 

It  appears  from  the  above  discussion  that  Grover's 
analytic  model  provides  an  improved  picture  of  the  quasi¬ 
static  KMP .  However,  the  fact  that  the  results  are  analytic 
places  some  restriction  on  the  answers.  Grover  points  out 
(Ref  2:21-24)  that  each  of  the  simplified  conductivity 
models  is  limited  by  time  and  range  values.  In  addition, 
the  polar  variation  of  the  source  terms  and  variations  in 
the  air  chemistry  parameters  are  not  considered. 


V. 


Late-Time  EMP  Sources  and 


Air  Chemistry  Parameters 

In  the  calculations  performed  by  Grover  and  the 
equivalent  numerical  results  presented  in  the  last  chapter, 
it  was  assumed  that  the  source  terms  (ionization  rate  and 
Compton  current)  wore  independent  of  the  polar  angle  theta 
and  the  air  chemistry  parameters  were  constant.  This  was  an 
obvious  requirement  to  obtain  analytic  results.  Monte  Carlo 
studies  have  shown,  however,  that  the  theta  dependence  of 
the  source  terms  can  be  significant,  especially,  for  ground 
capture  sources  (Refs  5  and  10) .  Furthermore,  the  electron 
air  chemistry  parameters  are  known  to  be  dependent  on  the 
total  electric  field  and,  also,  the  water  vapor  content  of 
the  air. 

In  this  chapter  the  late-time  EMP  sources  are  discussed. 
The  field  and  water  vapor  dependence  of  the  air  chemistry 
parameters,  also,  are  presented. 

Late-Time  Sources 

In  the  late-time  regime,  the  dominant  EMP  sources  result 
from  ground  and  air  capture  neutron  interactions  (Ref  6:13). 
The  gamma  rays  released  in  these  processes  interact  with  the 
atmosphere  producing  ionization  and  the  Compton  current  which 
drives  the  EMP. 

From  approximately  10-*4  to  5  x  10'3  seconds,  the  major 


gamma  source  conies  from  neutrons  which  are  captured  in  the 
ground  near  the  burst.  The  gammas  emitted  in  this  process 
are  angularly  dependent  due  to  the  presence  of  the  ground. 
Gammas  emitted  normal  to  the  surface  will  suffer  less 


attenuation  than  those  escaping  horizontally.  Consequently, 
the  ionization  rate  and  current  sources  resulting  from  this 
gamma  flux  will  vary  with  theta.  In  developing  curve  fits 
for  the  source  terms  from  Monte  Carlo  calculations,  O'Dell, 
et  al.  (Ref  10)  determined  that  the  ionization  rate  varied 
by  at  1  east  2.3  between  the  horizontal  and  vertical  direc¬ 
tions.  The  difference  was  found  to  be  even  greater  for 
larger  ranges.  More  significant  was  the  theta  dependence  of 
Jf.  Currents  pointing  straight  up  (0  =  0°)  were  as  much  as 
17  times  larger  than  the  currents  along  the  ground.  The 
variation  in  J  was  not  as  significant.  As  would  be  expected 
from  symmetry,  J„  approached  zero  for  0=0°  and  had  its 
maximum  value  at  0  =  90°.  To  model  this  a  theta  dependence 
of  (1  -  COSO)  was  used  (Ref  10:38). 

From  about  5  xl0~3  to  10" 1  seconds,  the  dominant  gamma 
ray  source  will  come  from  neutron  interactions  with  the 
atmosphere.  As  a  result,  the  ionization  rate  and  Compton 
current  are  nearly  independent  of  theta  (Ref  10:22).  O’Dell, 
et  al.  found  a  variation  of  1.3  in  the  ionization  rate  and 
2.3  for  the  radial  current.  Jn  was  found  to  be  much  smaller 
than  Jr  and  was  ignored. 

Shown  in  Figures  9,  10,  and  11  are  plots  of  the  ioniza¬ 
tion  rate,  the  radial  component  of  the  Compton  current,  and 
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Racial  Current  (Anys/m**2^ 


Radius  (m) 

Figure  10.  Radial  Compton  Current  (-J  )  vs  Radius  at  Four 
Angles,  Theta  Independent  Approximation  (Eq  2.11)  Shown 

for  Comparison 
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Figure  11.  Theta  Compton  Current  (J  )  vs 
Radius  at  Three  Angles  0 
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the  polar  component  (JQ)  at  several  different  angles.  The 
plots  were  generated  from  the  fits  described  above.  The 
curves  correspond  to  a  time  of  one  millisecond  after  a  10 
meqaton  surface  burs„.  For  comparison,  the  equations  used 
in  the  earlier  calculations  (see  Eqs  (2.11)  and  (2.22)]  are 
plotted  in  Fiqures  9  and  10.  It  can  be  seen  that  for  angles 
off  the  ground,  the  theta  independent  approximations  under¬ 
estimate  the  values  found  using  the  fits.  In  addition,  the 
theta  dependence  of  the  ionization  rate  is  a  function  of 
range  and,  as  shown  in  Figure  9,  varies  substantially  from 
the  approximate  result  for  distances  of  2000  meters  or  more. 
Note  that  falls  off  very  rapidly  with  distance  and  is  the 
largest  at  e  =  90°  as  expected. 

Figure  12  shows  the  time  dependence  of  J  and  J  at  a 

r  b 

particular  range  (1000  m) .  The  radial  component  is  given  at 
9  =  90°  and  0=0°.  As  discussed  earlier,  the  polar  varia¬ 
tion  of  J  reduces  with  time  as  the  dominate  source  switches 
r 

from  ground  to  air  capture  gammas.  JQ  is  seen  to  fall  off 
quite  rapidly  once  the  ground  capture  sources  are  gone.  The 
time  dependence  of  the  ionization  rate  is  the  same  as  that 
of  the  radial  current. 

The  curve  fits  outlined  above  were  developed  to  rep¬ 
resent  a  "typical"  thermonuclear  weapon.  They  were  biased 
toward  ranges  of  2  kilometers  or  less  partly  because  of 
statistical  uncertainties  in  the  Monte  Carlo  data  at  greater 
ranges.  Also,  this  is  the  range  of  interest  in  most  EMP 
calculations.  In  addition,  self-consistent  effects  between 
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the  fields  and  currents  were  not  included.  Although  this 
can  be  important,  the  approximation  is  probably  not  bad  at 
late  times  (Ref  10:36).  The  exact  form  of  the  fits  and  the 
numbers  used  to  qenorate  Figures  9-12  is  given  in 
Appendix  C. 

Air  Chemi s t ry  Parameters 

To  obtain  analytic  solutions,  the  air  chemistry 
parameters  used  to  determine  the  conductivity  were  assumed 
to  be  constant.  if  ions  dominate  the  conductivity ,  this  is 
a  good  approximation.  As  was  demonstrated  earlier,  however, 
electrons  outnumber  the  ions  at  shorter  ranges.  As  will  be 
seen  below,  the  electron  mobility  and  attachment  rate  are 
dependent  on  the  electric  field  and  water  vapor  content  of 
the  air.  Thus  for  points  close  to  the  burst,  the  calcula¬ 
tion  of  sigma  is  complicated  by  these  variations. 

This  problem  has  been  investigated  by  several  authors. 
In  general,  the  variations  in  i‘^  and  have  been  modeled 
by  curve  fits  to  the  available  experimental  data.  Many  EMP 
codes  employ  the  analytic  expressions  developed  by  Longley 
and  l.ongmire  (Ref  7).  Shown  in  Figure  13  is  the  fit  for  the 
electron  attachment  rate  as  a  function  of  the  total  field 
for  various  water  vapor  contents.  The  electron  mobility  is 
plotted  in  Figure  14.  The  curves  are  given  for  sea  level 
air  pressure  and  normal  density. 

Recently,  a  question  has  been  raised  as  to  the  accuracy 
of  the  data  from  which  Longley' s  fits  were  derived.  Based 
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Figure  11.  Attachment;  Rate  vs  Total  Electric  Field 

on  a  newer  set  of  experimental  data,  Pettus  and  Crevier 
(Kef  11)  developed  different  equations  for  the  mobility  and 
attachment  rate.  These  fits  are  slightly  more  complicated 
and  are  not  presented  here;  however,  Smith  and  Radasky  (Ref 
11)  showe'd  an  equation  which  demonstrated  the  effect  of  the 
new  data  on  the  electron  mobility.  The  equation  (Ref  13:30) 
limits  the  maximum  allowable  mobility  such  that 


1. 


max 


(1.0  +  96.77  (w)0,789) 


(m2/V-sec) 


(5.1) 


where  w  is  the  fraction  of  water  vapor  present.  This  modi- 
shown  in  Figure  14  for  the  one  percent  case. 
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Figure  14.  Mobility  vs  Total  Electric  Field 

It  is  obvious  that  these  two  parameters  are  not 
accurately  known.  The  debate  has  gone  on  for  a  number  of 
years  and  most  likely  will  continue.  Overall,  the  expres¬ 
sions  discussed  above  are  felt  to  be  reasonable  (Ref  3) . 
Moreover,  using  them  allows  variations  in  the  field  and  the 
water  vapor  content  to  be  considered.  The  exact  form  of 
bong  Icy  and  Longmire's  equations  is  given  in  Appendix  C. 
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VI .  An  Improved  Numerical  Method 


As  seen  in  the  previous  chapter,  the  ionization  rate 
and  Compton  current  sources  are  angularly  dependent  due  to 
the  presence'  of  the  ground.  Furthermore,  it  was  shown  that 
the  electron  mobility  and  attachment  rate  vary  with  the 
electric  field  and  the  water  vapor  content  of  the  air. 

These  v.ir  i  it  ions  con  Id  not  be  included  in  the  analytic 
expression?:  <>l  l.ongmiro,  Hill,  and  Wyatt  or  in  Grover's 
model.  A  loqical  extension  to  the  work  of  the  above- 
mentioned  authors  would  be  to  develop  a  solution  to  Maxwell's 
equations  which  could  incorporate  the  variations  in  sources 
and  air  chemistry  parameters. 

In  this  chapter  a  numerical  solution  will  be  presented 
which  allows  the  angular,  the  field,  and  the  water  vapor 
dependencies  to  be  included. 

S ol uti o n  Devel opmen t 

As  shown  earlier  (Chapter  II),  the  electric  field  in 
the  qua?;  i-s  tat  ic  phase  is  derivable  from  a  scalar  potential 
<f> .  Applying  this  definition  to  Maxwell's  equations  resulted 
in  the  following  differential  equation: 

-  —  ( r  2  E  ) - - - —  (SING  E  )  -  —  —  E 

r2  ;)r  rSINO  00  a  Or 

(6.1) 

-  —  —  E  =  —  —  ( r2 J  )  +  — - - —  (SIN0J  ) 

ar  00  or2  Or  rSINO  06 
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In  the  previous  solutions  to  this  equation,  and  o 
were  assumed  to  be  independent  of  theta  and  JQ  was  set  equal 
to  zero  (except  by  Wyatt).  Lonymire  and  Wyatt,  also, 
assumed  that  K  was  small  and  iynored  it.  Grover  included 
i:  bul  this  required  an  approximate  model  for  the  conduc¬ 
tivity.  These  restrictions  will  now  be  relaxed. 

Recalling  Eq  (3.1)  for  $(r,0) 


t(r,0)  -  Z  A  (r)P  (x) 
It  v  odd 


(6.2) 


Let  2  mean  2  .  Substituting  Eq  (6.2)  into  Eq  (6.1) 
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or  after  expansion  and  the  use  of  Eq  (3.3) 
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where  P^(x)  is  the  l'st  Associated  Legendre  polynomial  as 
before . 

Multiplying  both  sides  of  Eq  (6.4)  by  P  (x)  and 
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From  the  relationship  qiven  in  Eq  (3.6),  it  is  seen 
that  the  summation  in  the  first  term  of  Eq  (6.5)  reduces  to 
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Using  this  and  the  defining  Eqs  (6.6),  (6.7),  and 


(6.8),  Eq  (6.5)  becomes 
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Removing  the  terms  in  the  summations  where  m  =  l  and  moving 
the  rest  to  the  right-hand  side  of  the  equation  gives 
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At  this  [joint  all  derivatives  are  replaced  by  their 
finite  difference  approximations.  A  given  function  f (x)  is 
replaced  by 


f(x)  =  f(x  )  =  f 
n  n 


(6.12) 


where  n  is  the  n'th  point  where  f(x)  is  being  determined. 


i .  e  .  x 


n-Ax.  Likewise 


d  f (x)  fn+l  "  f n-1 


dx 


2Ax 


(6.13) 


47 


-  •  \<A. 

'W  - 


and 
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Applying  those  relationships  to  Eqs  (6.6)  and  (6.7) 


leads  to 
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P 

where  a  is  the  value  of  the  conductivity  at  the  n ’ th  radial 
n  1 

[joint  and  the  p  *  th  angle  theta. 

From  the  fits  for  JQ  (see  Appendix  C)  it  was  shown  that 


J0  =  C  •  (1  -  COS0) 


(6.17) 


where  C  is  a  function  of  radius  and  time;  hence 
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Using  this  result,  Eq  (6.8)  becomes 
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where  r^  =  n-Ar  and  is  the  radial  Compton  current  at  a 
given  range  point  n. 

Finally  Eq  (6.11)  becomes 
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Comparison  of  this  equation  with  Eq  (3.9)  reveals  that 
the  complex  derivation  above  has  resulted  in  an  equation 
whose  form  is  similar  to  the  earlier  result:  namely. 
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Consequently,  the  same  algorithm  that  was  presented  earlier 

can  be  used  to  determine  the  coefficients  An  . 

m 

The  solution,  however,  is  complicated  by  a  number  of 
factors.  In  the  previous  calculations,  the  coefficients 


a  ,  b 
m, n  m, n 


c  and  the  right-hand  term  R1?  could  easily 
m ,  n  d ,  m 

be  computed.  This  is  no  longer  the  case.  Since  the  conduc¬ 
tivity  and  Compton  current  are  allowed  to  vary  with  theta 
the  integrals  in  Eqs  (6.15),  (6.16),  and  (6.19)  are  much 

more  complicated  and  must  be  evaluated  numerically.  More¬ 
over,  the  term  R1?  includes  summations  which  involve  the 
'  d  ,m 
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Legendre  polynomial  coefficients  that  are  being  solved  for. 
Finally,  the  calculation  of  sigma  is  complicated  by  the 
field  dependence  of  the  electron  mobility  and  attachment 
rate . 

Solut ion  Technique 

To  solve  Kq  (6.22),  an  iterative  scheme  was  used.  An 
initial  guess  for  the  field  was  obtained  by  ignoring  the 
sununat Lon  terms  on  the  right-hand  side  of  the  equation  and 
using  average,  field  independent,  values  for  and 
Using  this  first  guess  as  the  input,  the  equation  was  solved 
again  but  the  summations  were  included  and  the  electron 
parameters  were  allowed  to  vary.  The  new  value  for  the 
total  field  obtained  here  was  then  compared  to  the  old  value 
at  all  range  [joints  and  at  one  angle.  This  scheme  was 
continued  until  the  difference  between  consecutive  itera¬ 
tions  was  less  than  a  specified  tolerance  at  all  the  range 
points.  It  was  assumed  that  convergence  along  one  angle 
guaranteed  convergence  everywhere. 

The  numerical  integrals  were  computed  using  Romberg 
integration  (Kef  1:206-211).  This  method  is  discussed  in 
Appendix  D . 

The1  resuLts  of  calculations  made  using  the  algorithm 
described  above  are  presented  in  the  next  chapter. 
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VII . 


Results 


Applying  the  algorithm  developed  in  the  last  chapter, 
the  quasi-static  electric  fields  were  computed.  The  ana¬ 
lytic  fits  for  the  source  terms  and  the  equations  for  the 
air  chemistry  parameters  presented  in  Chapter  V  were  used  in 
the  calculations.  The  results  are  presented  in  this  chapter 
along  with  a  comparison  to  the  earlier  calculations  given  in 
Chapter  IV.  The  effects  of  some  of  the  physical  parameters, 
also,  are  discussed. 

Test  Problem 

Usinq  the  same  yield  (10  Mt)  and  time  after  the  burst 
( 1 0 ~ 3  sec),  the  radial  and  polar  components  of  the  electric 
field  were  calculated.  The  radial  field,  Er,  is  plotted  as 
a  function  of  radius  in  Figure  15  and  polar  angle  in  Figure 
16.  It  is  seen  from  Figure  15  that  a  sizeable  radial  field 
was  found  for  points  close  to  the  burst.  Note  that  the 
field  becomes  negative  at  approximately  2000  meters.  This 
occurs  when  the  return  current  becomes  larger  than  the 
source  current.  The  polar  (theta)  field  is  given  in  Figures 
17  and  18.  A  peak  field  of  about  45  kV/m  was  obtained  at 
approximately  1200  meters  for  6  =  90°. 

The  difference  between  these  results  and  the  earlier 

calculations  was  significant.  Shown  in  Figures  15  and  17 

are  the  peak  radial,  E  (0°),  (dashed  line)  and  theta,  E 

r*  0 
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Figure  15.  Er  vs  Radius,  Full  Solution 
Initial  Solution  (Dashed  Line,  6=0°)  Also  Shown 
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Figure  16.  E  vs  Angle,  Full  Solution 
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e  17.  -  E0  vs  Radius,  Full  Solution.  Eq.  2.15,  (Dashed  Line 

and  Initial  Solution,  (Dot-Dashed  Line,  0  =  90°) ,  Also  Shown 
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(00°),  (dot  -dash  1  MU')  fields  that  wore-  computed  earlier 

us  mu  C.  row  i  1  s  sources  and  the  ,ur  chemistry  equations  for 

sigma  (s<-e  I  njures  4-7).  The  peak  fields  found  using  the 

improved  sol  u  t  ion  w.  ■  r  e  .  i[>pro.\  i  ma  t  e  1  y  1.5  times  larger  than 

the  pievnais  results.  In  fact,  the  maximum  K  found  with 

r 

the  new  solut  i on  was  almost  as  large  as  the  peak  from 
be  tore . 

The  polar  field  was  found  to  be  at  least  twice  as  large 
as  the  radial  !  ri'ld  at  their  peak  angles;  the  difference 
bein')  even  laruer  at  greater  distanees. 

Longmire’s  approx i matt*  solution  at  o  =  90°  is  shown 
for  comparison  in  Figure  17  (dish,  a  line).  Tt  is  seen  that 
the  overall  shape  of  the  tod,;  ;s  quite  different  than  the 
result,  obtained  in  the  calcul.it  ions  presented  here.  In 
addition,  it  should  be  noted  that  the  approximate  air 
chemistry  •  •nu.it  tons  were-  used  to  compute  s.iqma  in  Longmire's 
equation.  The  inherent  limitations  of  this  solution, 
therefore,  should  be  remembered. 


Theta  depend.  Tier  - 

It  was  found  in  all  calculations  that  the  dominant 
Legendre  polynomial  coefficient  was  Aj  (r).  Recalling  Eqs 
(3.17)  and  (  J.19) ,  one  would  expect,  therefore,  that 
should  behave  like  COSO  and  E  like  SIN0. 

From  Figure  16  it  can  be  seen  that  the  theta  dependence 
of  is  generally  like  COSO.  Note  that  at  larger  ranges 
becomes  nearly  theta  independent,  but  the  magnitude,  also, 
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is  much  smaller. 


F i <iure  18  shows,  however,  that  does  not  look  much 
like  S 1 N o .  Obviously  the  theta  dependence  is  complicated 
by  the  hOjher  order  polynomials.  The  theta  dependence 
predicted  by  I.onqmire's  solution  was  TAN(0/2),  see  Eq  2.15. 
This  equation  is  plotted  in  Figure  18  as  the  dashed  line  for 
r  “  1500  meters. 

it  was  found  in  the  calculations  that  five  polynomials 
were  sultieiont  to  properly  describe  the  theta  dependence. 
More  than  these  were  unnecessary  because  the  computed 
coefficients  were  much  smaller. 

boundary  Condition  at.  r  --  0 

To  find  the  values  of  the  integral  terms  in  Eqs  (6.15), 
(6.16)  ,  and  (6.19)  ,  it  was  necessary  to  know  J^.  and  o  at 
r  -•  0.  This  presented  a  problem  since  the  fits  used  to 
determine  these  parameters  had  a  1/r2  dependence.  Hence, 
at  r  -  0  J  and  o  would  be  infinite.  This  is  obviously  not 
possible  physically.  To  avoid  the  problem,  a  value  of 
r  -  6r  was  chosen  where  ir  was  allowed  to  approach  zero. 
Values  of  r  =  0.001  meters  or  smaller  resulted  in  the  same 
final  answers. 

Field  Dependence 

To  determine  the  effect  of  the  field  dependence  of  pe 
and  the  program  was  run  using  the  average  values 
•^resented  earlier  in  this  work,  i.e. 
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p  -  0.25  m'VV'sec 
e 

and 

a  =  1.5  x  108  sec- 1 
e 

The  results  are  shown  in  Figures  19  and  20.  For  comparison 

the  results  from  the  complete  solution  for  E  at  0°  and  E„ 

1  r  0 

at  90°  are  plotted  as  dashed  lines. 

The  variation  in  E  was  found  to  be  quite  large,  up  to 
25-30  percent.  The  difference  in  H  was  smaller,  but  as 
bag  as  B-LO  percent.  It  would  seem,  therefore,  that  the 
tio.ld  dependence  of  p  and  u  is  important. 

Effect  of  .1 

if 

In  t he  earlier  calculations  that  were  made  (Chapter  IV) 
no  theta  component  of  the  Compton  current  was  included.  It 
was  shown  that  J(|  was  much  smaller  than  J^.  except  for  r  <  \ 
(Kef  10).  To  determine  its  effect,  a  calculation  was  made 
where  J  was  set  equal  to  zero.  The  results  are  shown  in 
Figures  21  and  22.  There  was  no  difference  in  the  radial 
field,  and  as  can  be  seen  from  Figure  22,  the  variation  in 
E  was  very  small.  The  difference  was  loss  than  5  percent 
at  500  meters  and  the  answers  were  the  same  for  ranges 
beyond  approximately  1500  meters. 

Variation  in  Electron  Mobi 1 i ty 

As  was  discussed  in  Chapter  V,  recent  data  has  shown 
that  the  fit  for  mobility  and  attachment  rate  developed  by 
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Figure  19.  Er  vs  Radius,  Field  Independent 
Full  Solution  Shown  as  a  Dashed  Line  (6=0° 
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Figure  20.  -  E0  vs  Radius,  Field  Independent 

Full  Solution  Shown  as  a  Dashed  Line  (0  =  90° 
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Figure  22.  -  E0  vs  Radius,  JQ  = 

Full  Solution  Shown  as  a  Dashed  Line  (e 


Long ley  and  Longmire  might  not  be  correct  (Ref  11) .  In  the 
calculations  made  in  this  report,  the  Longley  fits  were  used 
along  with  the  modification  to  the  mobility  suggested  by 
Smith  (Ref  13:30).  To  see  what  effect  this  modification  had 
on  the  fields,  a  calculation  was  made  which  did  not  include 
Smith's  equation.  The  results  are  given  in  Figures  23  and 
24.  Once  again  the  original  solutions  for  E^  at  8  =  0° 
and  Et.  at  0  =  90°  are  shown  for  comparison  (dashed  lines)  . 
The  variation  in  and  EQ  was  less  than  2  percent  for  the 
peak  field  values.  Obviously  the  modification  makes  some 
difference;  however,  the  air  chemistry  parameters  are  not 
known  with  enough  certainty  to  suggest  any  major  conclusion 
(Ref  3:29).  It  should  be  pointed  out,  though,  that  no 
change  was  made  in  the  electron  attachment  rate  equations. 

Summary 

The  results  of  this  section  suggest  that  the  polar 
variation  of  the  ionization  rate  and  Compton  current  are 
important  for  the  calculation  of  the  electric  fields.  In 
addition,  it  was  seen  that  variations  in  the  air  chemistry 
parameters  and  the  theta  component  of  the  Compton  current 
have  some  effect  on  the  solutions.  In  the  next  chapter  a 
series  of  parametric  studies  will  be  presented  which  will 
further  show  the  importance  of  these  parameters. 
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Figure  23.  E  vs  Radius,  Electron  Mobility  Not  Limited 
Full  Solution  Shown  as  a  Dashed  Line  (6  =  0°) 
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VIII  . 


Parametric  Studies 


Usirtq  the  code  developed  from  the  algorithm  in  Chapter 
VI,  the  impact  of  changing  various  parameters  was  consid¬ 
ered.  Both  the  accuracy  of  the  code  and  the  effects  of 
varying  physical  parameters  were  studied.  The  results  are 
presented  in  this  chapter. 

Code  Accuracy 

The  computer  program  used  in  this  report  (see  Appendix 
F)  was  tested  by  varying  the  numerical  parameters  which 
affected  the  calculations. 

In  approximating  the  derivatives  by  finite  differences, 
a  Ar  had  to  be  chosen.  It  was  found  that  for  a  total  range 
of  0-5000  meters  a  Ar  of  50  meters  was  sufficient  to 
guarantee  convergence  in  the  solutions;  hence  only  100  range 
points  were  needed.  Of  course,  more  could  have  been  used 
but  the  added  memory  requirements  and  calculation  time  were 
unnecessary .  A  AO  of  one  degree  was  used  for  the  derivative 
calculations  with  respect  to  o . 

As  discussed  in  Appendix  l),  the  Romberg  integration 
technique  does  not  use  a  predetermined  step  size  when  the 
integrals  are  calculated.  Rather,  a  so-called  "row  limit" 
is  selected  which  allows  for  any  desired  accuracy.  It  was 
found  that  a  row  limit  of  four  was  sufficient  for  calcu¬ 
lation  of  the  numerical  integrals.  Choosing  five  gave 


sliqhtly  bettor  accuracy,  but  the  calculation  time  for  a 
complete  solution  was  doubled. 

In  the  iteration  scheme  used  to  find  the  fields,  a 
tolerance  of  100  V/m  was  selected  for  the  convergence  test 
between  two  consecutive  calculations .  In  addition,  it  was 
determined  that  the  answers  were  the  same  regardless  of  the 
angle  selected  for  the  convergence  test. 

The  final  numerical  parameter  considered  was  the  number 
of  Legendre  polynomial  coefficients,  A  (r),  that  were 
necessary  to  adequately  describe  the  potential  function 
<t>(r,i3)  •  It  was  determined  that  the  first  five  odd  terms  in 
the  series  were  sufficient  (i.e.  i  =  ],  3,  5,  7,  &  9). 

Using  six  terms  did  not  improve  the  answer  and  required 
twice  as  much  computer  time  for  a  full  calculation. 

'The  parameters  discussed  above  are  summarized  in  Table 
I  below. 


Table  1 

Computer  Code  Parameters  Used  in  Calculations 


Number  ol  Range  Points  -  100 

A  r  -  50m 

A0  -  1  deg 

Romberg  Row  L imit  -  4 

Convergence  Tolerance  -  100  V/m 

Number  of  Polynomials  -  1st  5  odd 


Tunc  Dependence 

Shown  in  Figure  25  is  a  plot  of  the  total  field  along 
the  ground  (0  --  90°)  vs  time  for  five  different  ranges. 


68 


LEGEND 
5000  meters 


Time  (sec) 

Total  Field  vs  Time 


This  field  corresponds  to  the  maximum  theta  component  (i.e. 

I-:  =  0).  The  two  plateaus  appear  to  agree  quite  well  with 

the  time  regimes  where  either  ground  or  air  capture  sources 
dominate.  The  magnitude  of  the  field  is  seen  to  decrease 
approximately  by  a  factor  of  ten  over  the  time  range  shown. 
This  is  an  interesting  result  since  the  current  was  shown  to 
fall  off  by  two  orders  of  magnitude  for  the  same  time  range 
(see  figure  12). 

figure  2f  shows  the  total  field  at  0  =  0°  (E  =  0)  . 

0 

The  time  dependence  is  approximately  the  same  as  above, 

although,  the  decrease  in  the  total  field  is  less  than  an 

order  i»l  magnitude.  If  is  seen  that  lor  times  greater  than 

approximately  10--'  seconds  the'  peak  radial  field  at  500 

meters  is  nearly  as  large  as  the  maximum  polar  field.  This 

result  differs  substantially  from  a  curve  given  in  a  report 

by  l.ongmire  (Ref  9:45)  .  The  calculations  presented  in  this 

report  show  a  much  larger  radial  field  than  shown  in  the 

above- re lerenced  curve  for  E  . 

r 

Additional  plots  for  the  total  field  at  10-2  and  10-1 
seconds  are  given  in  Appendix  E. 

Y  i  e  U1  Depend*  Mice 

figures  27  and  28  show  the  yield  dependence  of  the 
total  Meld  at  a  90°  and  "  0"  respectively.  In 

figure  27  it  is  seen  that  the  total  field  is  nearly  in¬ 
dependent  of  the  yield  at  500  meters  and  falls  off  slightly 
for  lower  yields.  for  larger  distances  the  total  field  is 
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Figure  26.  Total  Field  vs  Time, 
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Figure  27.  Total  Field  vs  Yield 
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Figure  28.  Total  Field  vs  Yield 


more  dependent  on  the  yield.  This  can  be  explained  by 
considering  the  ratio  of  J/o.  I-'or  short  ranges  this  ratio 
is  nearLy  independent  of  yield  because  electrons  dominate 
the  conduct i vi ty .  At  greater  distances  the  conductivity  is 
dominated  by  ions  and  the  ratio  is  no  longer  yield  inde¬ 
pendent.  A  similar  yield  dependence  is  shown  in  Figure  28 
for  the  total  1 ield  at  o  -  0°. 

Additional  curves  for  1000  KT  and  100  KT  are  given  in 
Appendix  K. 

Water  Vapor  Dependence 

The  linal  parameter  considered  was  the  effect  of 
various  water  vapor  contents  on  the  total  field.  In  Figures 
29  and  30  the  total  field  is  plotted  as  a  function  of  water 
vapor  content  in  the  air.  As  seen  in  Figure  29,  the  field 
at  0  -  90°  varies  quite  substantially  over  the  range  of 
water  vapor  fractions  shown.  Recalling  Figures  13  and  14 
for  the  mobility  and  attachment  rate,  it  was  seen  that 
higher  water  vapor  contents  decreased  the  mobility  and  in¬ 
creased  the  attachment  rate.  Thus,  as  the  fraction  of  water 
vapor  in  the  air  increases  the  conductivity  decreases. 

Since  E  -  J/o,  the  field  should  get  larger.  This  result  is 
seen  clearly  in  Figure  29.  Note  that  the  dependence  is 
greater  for  the  500  meter  curve  due  to  electron  dominated 
conductivity.  For  larger  ranges  where  ions  dominate,  the 
dependence  is  slightly  less.  Figure  30  shows  that  the  field 
at  0  -  0°  does  become  nearly  independent  of  the  water  vapor 
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Figure  29.  Total  Field  vs  Water  Vapor  Content,  6  =  90° 


content  for  ranges  of  1500  meters  or  larger. 

Additional  curves  for  water  vapor  fractions  of  0.00 
and  0.04  are  given  in  Appendix  E. 


IX.  Conclusions  and  Recommendations 

Cone  1  us i uns 

In  this  study  the  quasi-statie  EMP  fields  were 
calculated  using  a  time-independent  numerical  method.  On 
the  basis  of  the  research  and  the  calculations  performed, 
the  following  conclusions  are  made: 

1 .  The  angular  dependence  of  the  ionization 
rate  and  Compton  current  are  important  in  the 
calculation  of  the  fields.  The  peak  field  values 
calculated  in  this  work  varied  by  as  much  as  50 
percent  from  Grover's  theta  independent  calcula¬ 
tions  (see  Figures  15-18).  Thus,  the  results 
presented  in  this  paper  extend  Grover's  work  by 
using  more  realistic  sources. 

2.  Even  though  the  radial  field  is  small  at 
late  times,  its  effect  is  important  for  a  complete 
solution  to  Maxwell's  equations.  In  Longmire's 
analytic  development  (Ref  9)  E  was  ignored  and 
the  resulting  EQ  field  had  a  limited  range  of 
validity.  lsy  including  H  ,  the  calculations  made 
in  this  work  provide  an  improved  picture  of  the 
late  time  EMP. 

3.  The  electric  fields  that  result  from  a 
surface  nuclear  burst  are  dependent  on  variations 
in  the  air  chemistry  parameters.  Using  field 
independent,  average  values  for  the  electron 
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mobility  and  attachment  rate  resulted  in  peak 
field  values  which  were  10  to  30  percent  different 
from  the  calculations  where  these  parameters  were 
allowed  to  vary.  As  shown  in  Figure  29  the  field 
is  also  a  strong  function  of  the  water  vapor 
content  of  the  air. 

4.  The  inclusion  of  the  theta  component  of 
the  Compton  current  (J  )  is  not  important  for  EMI' 
calculations  as  it  only  changes  the  fields  for 
points  close  to  the  burst.  In  the  test  case 
studied  in  this  work,  the  variation  was  less  than 
5  percent  for  ranges  shorter  than  500  meters.  The 
answers  were  identical  for  larger  ranges. 

Recommendations 

Based  on  the  results  obtained  in  this  study,  the 
following  recommendations  are  suggested: 

1.  The  results  of  the  time-independent 
calculations  presented  in  this  report  should  be 
compared  to  the  results  obtained  with  time- 
dependent  codes  (such  as  LEMP)  to  determine  how 
well  the  fields  are  being  predicted. 

2.  The  code  should  be  modified  to  include 
the  new  fits  (Ref  11)  for  the  electron  mobility 
and  attachment  rate.  A  comparison  should  then  be 
made  with  the  results  presented  in  this  work. 

3.  The  validity  of  the  quasi-static 


approximation  should  be  considered. 
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Appendix  A 

I '  mite  DH'Icrmh'i's  in  Initial  Sol  ut  ion 

In  thin  appendix  the  tinite  difference  approximations 
which  Led  to  Kii  (1.9)  .ire  presented. 

From  Maxwell’s  equations  and  the  Ljeqendre  polynomial 
expansion  the  following  difierential  equation,  Fq  (3.8), 
was  ob t  .  i  i  in 'd  : 


i  •  [  i  ; 

;i  t 


I  .in 

o  :.i  r 


A  (r) 
l(t+l)  -  — — 
r  2 


CMt+l) 


a  i 


( x )  d  x 


(A-l) 


To  obtain  a  numerical  solution,  all  the  derivatives  in 
llq  (A-l)  are  replaced  by  their  central  difference  approxi¬ 
mations.  Lot  r  -  n*Ar  where  n  is  an  integer.  Hence,  for 
n 

N  range  points,  Ar  is  found  from  Ar  =  rmlx/N .  Then  for  a 
given  t  unit  ion  f(r),  the  following  are  obtained: 

f(r)  -  f(rn)  -  fn  (A-2) 


A 1  so 


d  r  (r) 
dr 


f  .  -  f  , 
n  +  1 _ n- 1 

2Ar 


(A-3) 


and 


d  ’  f  (  r )  __  n-tl  ^  ^n  +  ^n-  1 

dx;>  Ar2 


( A-4 ) 
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This  appendix  contains  the  results  of  calculations  that 
were  made  using  Grover's  sources  and  his  Model  (3)  conduc¬ 
tivity.  These  results  are  discussed  in  Chapter  IV.  The 
figures  are  grouped  as  follows: 

Figures  31  -  34  ;  Fields  found  using  oq  =  178.0 

Figures  35  -  38  ;  Fields  found  using  oq  =  143.2 

Note  that  in  all  figures  the  absolute  value  of  the 
field  is  plotted.  Actually  is  positive  and  E0  is 
negative  over  most  of  the  range. 
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Figure  31.  E  vs  Radius,  Simplified  Sources 
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Figure  32.  E  vs  Angle,  Simplified  Sources  oQ  =  178.0 
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Figure  33.  -  E0  vs  Radius,  Simplified  Sources 
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Figure  34.  -  E0  vs  Angle,  Simplified  Sources 
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Figure  36.  E  vs  Angle,  Simplified  Sources  <tq  =  143. 
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Appendix  C 

Analytic  Fits  for  Sources  and 
Ai r  Chemistry  Parameters 

In  this  appendix  the  analytic  equations  for  the  source 
terms  and  air  chemistry  parameters  are  presented.  These 
fits  were  discussed  in  Chapter  V  and  were  used  in  the 
calculations  presented  in  Chapters  VII  and  VIII. 

Source  'terms 

The  source  functions  are  given  in  Table  TI  (Ref  10:38) 
The  symbols  and  units  used  in  these  formulas  are  as  follows 

r  =  radial  distance  from  the  burst  (cm) 
p  =  COSO  ;  6  =  polar  angle 

Y  =  total  yield  (KT) 

NtNa  =  number  of  neutrons  per  kiloton 
p  =  air  density  (mg/cm3) 
t  =  time  (seconds) 

Y  =  ionization  rate  (ion  pairs/cm3 • sec) 

Jr  =  radial  current  density  (abamps/cm2 ) 

-  polar  (theta)  current  density  (abamps/cm2) 

In  all  calculations  the  following  values  were  assumed: 

N  N  =  2.0  x  10 2 3  neutrons/KT 

t  cl 

p  =  1.225  mg/cm3 


Table  II 


<<  J  for  all  r  and  y  of  interest 


Note  also  that  y,  and  were  converted  to  MRS  units  for 
use  in  the  program. 

Air  Chemistry  Parameters 

To  determine  the  electron  mobility  and  attachment  rate, 
the  analytic  equations  developed  by  Lonqley  and  Longmire 
(Ref  7)  were  used  except  for  one  modification.  The  maximum 
mobility  for  a  given  water  vapor  content  was  found  using 
Smith's  equation,  Kq  (5.1).  Jn  all  cases  the  calculations 
were  made  using  whatever  units  the  fits  were  in  and  at  the 
end  the  answers  were  converted  to  MRS  units.  The  equations 
are  summarized  in  Table  III.  The  symbols  and  units  used  in 
these  formulas  are  as  follows: 

l  =  E/p  (electrostatic  units/atmosphere) 
r q  =  esu/atm  ;  found  from  cQ  =  f(e,w) 
w  =  fraction  of  water  vapor  present 
P  -  percent  water  vapor  (P  =  100 -w) 

=  electron  mobility  (cm/sec*esu) 

Pr  =  p  /pQ  =  ratio  of  actual  air  density  to  normal 
=  three-body  attachment  rate  (sec-1) 

** 2  ^  ~  two-body  attachment  rate  (sec-1) 

In  all  calculations  the  following  were  assumed: 

pq  =  1.255  x  10-3  gm/cm3 
1  esu  =  30,020.0  Volts/m 
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Table  III 


Mobility  and  Attachment  Rate  Equations 


The  Mobi 1 i ty 

t  i  E/p  (esu/atmosphere)  P  s  100  w  k 


* U  '  1  t  Ak  '  ,or  e 


r)' 

-  C  k  ,  (or  l  >  l  . 


+  t 


Bk 

2 


A 

B 

C 


for  c j  < e  < r 


0.834 

P 

2.457 

0.6884 

1.195 

0.07853(1  +  Ak) 
3.015  +  Ck 


/  cm  v  _  1 . 0  x  1 0 6 

’  16.8  +  cQ 

\ sec  esu  /  P 

0.6  3  +  26 . 7 

r 

o 

0.60 


Pa(V 

We ^  C  ^  1  -  w  +  wR ( eg ) 


R  (  Gg  )  =  1.55  + 


210 


1 +  11.8e0  + 7 . 2gq 


Ttie  A  t  tac’htiK-nt  Rate 
throe-body 


a,  (sec-1)  =  108p2 
3a  r 


0.62  +  800e2 


1  +  103e2[Co(i  +  0 . 03e2) J 1/3 


two-body 


(sec-1)  =  1.22  x  108p  e  -2 1 • 15/e( 
2a  r 


<*e  =  (1  -  w){(l  +  34 . 4w)  a^a  +  «>2a^ 
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Appendix  D 
Romberg  7 ntegrat ion 


To  find  the  numerical  integrals  in  the  improved 
solution,  Romberg  integration  was  used.  This  method  is 
briefly  outlined  below. 

Romberg  integration  starts  by  using  the  trapezoid  rule 
to  give  a  first  approximation  to  the  integral.  The  general 
formula  for  trapezoidal  integration  is 


b 

J  f  ( x )  d  X 

a 


h 


f (a)  +  1(b) 


+ 


n-1 
2  i 
J-l 


1 


f  (x^ 


+  0(h2) 


(D-l) 


where  h  b-a/n  and  x^  =  a  +  j • h  for  each  j  =  0,1,..., n. 

The  error  term  is  0(h2).  By  applying  Richardson  extrapo¬ 
lation  (Ref  1:208)  an  improved  expression  is  obtained  which 
converges  faster  and  has  an  error  term  of  0(h4)  .  The 
equation  is 


/ 


f (x) dx 


4R  .  -  R.  .  . 

— Lt 1 - Jizlil  +  o  (h4  ) 


(D-2) 


The  R,  ,  terms  are  found  from 
k ,  1 

hl 

If  (a)  +  f  (b)  J 

i ,  1  2 


(D-3) 


where  h^  =  b-a,  and  in  general 


i 


Rk-l,l  +  hk-l  1  f(a+  ll'2>  hk-l> 

(D-4 ) 

l  i  =  1  J 

in  which  =  j  • 


Finally,  it  can  be  shown  that  the  technique  is 
generalized  to  any  R^  ^  such  that 


J  f(x)  *  R 


where 


R  . 
1 


/  3 


4^-1  R. 


i/  3-1 


i  3*1 


-  R.  -  .  , 

i-l, 3-1 

1 


(D- 


for  eacli  i  =■  2,3,4,...,n  and  j  =  l,2,...,i.  The  approxima 
tions  R ^  .  can  be  represented  as  in  Figure  39. 


Rl,l 

R2  , 1 

R2.2 

H3,l 

R  3 , 2 

R3,3 

H4,l 

R4,2 

R4.3 

R4,4 

• 

• 

• 

• 

• 

• 

Rn,l 

Rn ,  2 

R  , 
n,  3 

Rn ,  4 

• 

• 

• 

R 

n,n 

Figure  39.  Generalized  Romberg  Terms 


Hence  the  integral  is  seen  to  converge  diagonally. 
Note  that  no  predetermined  step  size  (h)  is  chosen.  Thus 


the  integral  is  computed  by  selecting  a  maximum  "row  limit" 

n  and  determining  all  the  terms  in  Figure  39  until  R  is 

n,n 

found.  The  method  is  fast  because  once  R.  1  is  known,  all 

1  y  1 

other  terms  in  the  row  are  found  using  this  and  the  value  of 
from  the  previous  row.  The  complete  algorithm  which 
was  built  into  the  computer  program  is  given  in  Figure  40. 


Romberg  Algorithm 

Select  a  positive  integer  n.  Let  h.  =  b-a 
and  perform  the  following: 


1. 

Set 

i  -  1 

2. 

Add 

1  to  i 

3. 

I  f  i 

>  n 

4  . 

Set 

R  •  , 

i,  1 

1 

2 

5. 

Add 

1  to  j 

f> . 

1  f  ) 

-  i  , 

7  . 

Set 

1,1  2 


[f (a)  +  f (b)  ] 


2  1 

+  h.  .  E  f  (a  +  (k-^)h  . ) 

1-1  k=l  ^  1“A 


R.  . 

i  .3 


2 2  (  i  1  *  R .  .  .  -  R .  .  . 

i>3-l  l-l. 3~1 


22  (i-D  _  x 


8.  C.o  to  Step  5 


9.  The  procedure  is  complete  and  R^  n  approximates 


f  f (x) dx 
a 


Figure  40.  Romberg  Algorithm  (Ref  1:210) 
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Appendix  E 
Parametric  Studies 


This  appendix  contains  the  results  of  calculations  made 
using  the  complete  solution  algorithm.  The  plots  represent 
a  series  of  parametric  studies  where  the  time,  the  yield, 
and  the  water  vapor  content  were  varied.  The  figures  are 
grouped  according  to  the  parameter  studied. 


Time  Dependence  (t) 


Figures 

41  - 

42  ; 

Data 

for 

t  = 

10 

2  sec 

Figures 

43  - 

44  ; 

Data 

for 

t  = 

10 

- 1  sec 

Yield  Dependence 

(Yld) 

Figures 

45  - 

46  ; 

Data 

for 

Yld 

= 

1000  KT 

Figures 

47  - 

48  ; 

Data 

for 

Yld 

= 

100  KT 

Water  Vapor  Content  (Wvc) 

Figures 

49  - 

o 

in 

Data 

for 

Wvc 

= 

0.00 

Figures 

51  - 

52  ; 

Data 

for 

Wvc 

= 

0.04 

The  total  field  was  found  by  taking  the  square  root  of 


E  2  +  E  2 
r  u 


( 
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Total  Field  vs  Radius,  Full  Solution  t  =  10“2  sec 
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Figure  42.  Total  Field  vs  Angle,  Full  Solution  t  =  10-2  sec 
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Figure  43.  Total  Field  vs  Radius,  Full  Solution  t  =  10"1  sec 
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Figure  44.  Total  Field  vs  Angle,  Full  Solution  t  =  10'1  sec 
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Figure  45.  Total  Field  vs  Radius,  Full  Solution  Yld  =  1000  KT 
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Figure  46.  Total  Field  vs  Angle,  Full  Solution  Yld  =  1000  KT 
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Figure  48.  Total  Field  vs  Angle,  Full  Solution  Yld  =  100  KT 
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Figure  49.  Total  Field  vs  Radius,  Full  Solution  Wvc  =  0.00 
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Figure  50.  Total  Field  vs  Angle,  Full  Solut 
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Figure  51.  Total  Field  vs  Radius,  Full  Solution  Wvc  =  0.04 
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Figure  52.  Total  Field  vs  Angle,  Full  Solution  Wvc  =  0.04 
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Program  Poser i ption  and  Li  sting 

The  program  used  to  find  the  late-time  EMP  was  written 
in  FORTRAN  V.  Every  effort  was  made  to  include  descriptive 
comments  and  documentation.  Structured  programming  tech¬ 
niques  wore  applied  and  many  of  the  calculations  are  made  in 
various  subroutines.  The  program  complies  with  the  1977 
ANSI  standard  and  should,  therefore,  be  transportable. 

Other  than  possible  changes  in  the  numerical  accuracy 
of  the  calculations  the  only  input  to  the  program  is  the 
time,  the  yield,  the  water  vapor  content,  and  the  relative 
air  density.  The  atmospheric  pressure  was  taken  to  be  one 
atmosphere  (sea-level) . 

A  flow  diagram  for  the  program  is  given  in  Figure  53. 
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Figure  53.  Program  Flow  Diagram 
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